COHERENT STRUCTURES AND CARRIER SHOCKS IN THE NONLINEAR 

PERIODIC MAXWELL EQUATIONS 



GIDEON SIMPSON AND M.I. WEINSTEIN 

Abstract. We consider the one-dimensional propagation of electromagnetic waves in a weakly nonlinear 
and low-contrast spatially inhomogeneous medium with no energy dissipation. We focus on the case of 
a periodic medium, in which dispersion enters only through the (Floquet-Bloch) spectral band dispersion 
associated with the periodic structure; chromatic dispersion ( time-nonlocality of the polarization) is ne- 
glected. Numerical simulations show that for initial conditions of wave-packet type (a plane wave of fixed 
carrier frequency multiplied by a slow varying, spatially localized function) very long-lived spatially localized 
coherent soliton-like structures emerge, whose character is that of a slowly varying envelope of a train of 
shocks. We call this structure an envelope carrier-shock train. 

The structure of the solution violates the oft-assumed nearly monochromatic wave packet structure, 
whose envelope is governed by the nonlinear coupled mode equations (NLCME). The inconsistency and 
inaccuracy of NLCME lies in the neglect of all (infinitely many) resonances except for the principle resonance 
induced by the initial carrier frequency. We derive, via a nonlinear geometrical optics expansion, a system of 
nonlocal integro-diffcrcntial equations governing the coupled evolution of backward and forward propagating 
waves. These equations incorporate effects of all resonances. In a periodic medium, these equations may 
be expressed as a system of infinitely many coupled mode equations, which we call the extended nonlinear 
coupled mode system (xNLCME). Truncating xNLCME to include only the principle resonances leads to 
the classical NLCME. 

Numerical simulations of xNLCME demonstrate that it captures both large scale features, related to third 
harmonic generation, and fine scale carrier shocks features of the nonlinear periodic Maxwell equations. 
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1. Overview 



Realized and potential applications of microstructured dielectric media motivate a thorough mathemat- 
cal study of wave-propagation governed by nonlinear hyperbolic equations, e.g. Maxwell's equations with 



periodic and nonlinear constituitive laws. This paper explores a class of nonlinear hyperbolic equations with 
a spatially periodic flux function: 



(1.1a) d t v + dj(x,v) = 

(1.1b) f(z,0)=0, 

(1.1c) f(l + 27T,v) =f(l,v). 

In particular, we shall assume that periodic variations are weak (a low contrast structure) and study solutions, 
whose amplitude is small and such that the effects of periodicity-induced dispersion and nonlinearity are in 
balance. 



Indeed, a non-trivial spatially periodic structure is dispersive. This can be seen by linearizing ( 1.1 1 about 
the state v = 0, giving the linear system: 

(1.2) d t V + d x (D v {(x,0)V)=0, 



which retains periodicity. Floquet-Bloch theory [7,38 implies that associated to the PDE (1.2) is a family of 



band dispersion functions k h-> Wj (k) , fee (— §, ^JT Wave propagation is dispersive since the group velocities, 
uj'j(k), are typically non-zero. Thus, waves of different wavelengths travel with different speeds. 

Dispersive properties, encoded in the functions Wj(-) and the associated Floquet-Bloch states, can be 
manipulated by tuning the periodic structure through, for example, modification of the periodic lattice, the 
maximum and minimum variations of D v f(x, 0) (material contrast), etc. 

It is well-known that for general initial conditions, solutions of hyperbolic systems of conservation laws 
with spatially homogeneous nonlinear fluxes: 

(1.3) d t v +dj(v) = 



develop singularities (shocks) in finite time, 20 22 



Question 1: Is spectral band dispersion, due to a periodic structure, sufficient to arrest shock formation?^ 

The ability to control or inhibit the formation of singularities in nonlinear wave propagation could have 
significant impact in, for example, electromagnetics and elasticity. Strictly speaking, the answer to Question 



1 is no. Indeed, for a system of the form (1.1), let us suppose that the flux function was periodically 
piecewise constant. Finite propagation speed considerations imply that for appropriate initial data, which 
are sufficiently localized within a uniform region, a shock will form. The dispersive character of the periodic 
structure is manifested only on sufficiently large spatial and temporal scales. Thus, the problem of controlling 
shock formation should be posed relative to some class of initial conditions. 

A second motivation is the study and design of media which support the propagation of stable soliton-likc 
pulses. These have applications to optical devices which transfer store or, in general, process information 
which is encoded as light pulses. Associated with dispersive wave-propagation at wavenumber fc* is a dis- 
persion length ~ (^"(fefr)) - 1 . Soliton formation is possible on length scales where the dispersion length and 
the characteristic length on which nonlinear effects act are comparable. Technological advances have made 
it possible to fabricate microstructured media with specified dispersion lengths at specified wavelengths. For 
a given dispersion length, a balance between dispersion and nonlinearity is achieved by tuning the strength 
of the nonlinear effects through adjusting the field intensity (by an amount which is material dependent). 
An example of this balance at work is gap soliton formation in periodic structures. These are experiments in 
optical fiber periodic structures (gratings) involving highly intense (nonlinear) light with carrier wave-length 
satisfying the Bragg (resonance) condition. The length-scale of such solitons is 10~ 2 meters [§]. 

Theory predicts the existence of gap solitons traveling at any speed, v, between zero and the speed of light, 
c [l][3]. Experiments [8] demonstrate speeds as low as .3c to .5c. Potential applications of gap solitons, based 
on the design of appropriate localized defects in a periodic structure, are all-optical storage devices [10]. 



^For example, though typical smooth initial data for the inviscid Burgers equation dtu + ud x u = develop shocks in finite 
time, the corresponding solutions of the Korteweg - de Vries (KdV) equation, dtu + ud x u + d x u = 0, a dispersive perturbation, 
remain smooth for all time. 
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The term gap soliton is used due to frequency of the gap soliton envelope lying in the spectral gap of the 
linearized system. 



Physical predictions of gap solitons are based on explicit solutions of nonlinear coupled mode equations 
(NLCME), given below in ( |2.9[ ). NLCME has been formally derived in, for example, [I] from (1.4); see also 
the discussion in Section [2j Rigorous derivations of NLCME, from models with appropriate dispersion have 
been presented for the anharmonic Maxwell-Lorentz equations 
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and other nonlinear dispersive equations; 



see e.g. 

Within the approximation of a small amplitude wave field as a wave-packet with slowly varying envelope 
and single carrier frequency, propagating through a low contrast periodic structure near the Bragg resonance 
( see scaling in Figure [lj ) , NLCME is argued to govern the principle forward and backward slowly varying 
envelopes of carrier waves; see [I] and references therein. 

As discussed in fllf and in section [1| if the only source of dispersion is the spatial dispersion of the 
periodic medium (e.g. negligible chromatic dispersion) for weakly nonlinear waves in low contrast media all 
nonlinearity- generated harmonics are resonant and therefore all mode amplitudes are coupled at leading order. 
The correct mathematical description would appear to require infinitely many interacting modes. Thus, the 
classical NLCME are not a mathematically consistent approximation. NLCME may however be satisfactory 
physical description, for some purposes. t Indeed, the soliton wave form prediction based on NLCME appears 
to describe some features of experiment. ^\ 



Question 2: Do nonlinear periodic hyperbolic systems have stable coherent structures, and can one develop 
a mathematical theory? How are the classical NLCME related to this theory? See the discussion in section 

m 



Physicists argue in two ways that the coupling to higher harmonics is argued to be negligible : (i) The material systems 
considered are dissipative at higher wave numbers. Higher wave numbers are damped and therefore these mode amplitudes can 
be ignored, and (ii) Chromatic dispersion (arising due to the finite time response of the medium to the field) causes nonlinearly 
generated harmonics to be off-resonance. Therefore, an initial condition exciting the principle modes will not appreciable excite 
higher harmonics. These rationales are somewhat ad hoc since the precise damping mechanisms are not well-understood and 
chromatic dispersion is a much weaker effect than photonic band dispersion for weak periodic structures. 
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(a) (b) 

Figure 2. On the left is a simulation of the Maxwell equations. On the right is the 
simulation of a truncated asymptotic system, resolving the first and third harmonics. Both 
simulations were initiated with the same initial conditions. The two side pulses about the 
main wave appear to be the result of third harmonic generation. 



In this article we report on progress on Questions 1 and 2 in the context of the one-dimensional, nonlinear 
Maxwell equations governing the electric (E) and magnetic (B) fields: 

(1.4a) 8 t D = d z B, 

(1.4b) d t B = d z E. 

with constitutive law 

(1.5) D = e(z,E) E = (n 2 (z) + X E 2 ) E 

(1.6) n(z) = n + eN(z) 

n(z) is a linear refractive index, consisting of a nonzero background average part, no, and a fluctuating (e.g. 
periodic) part eN(z). The nonlinear term \E 2 is the nonlinear refractive index, arising from the Kerr effect; 



in regions of high intensity the refractive index is higher. The consituitive law (1.5), prescribes D as a a 
local function of E. Thus chromatic dispersion, which arises due a time-nonlocal relation between D and E 
has been neglected. For simplicity, we assume ng = 1, which can be arranged by a simple scaling. 

1.1. Summary of results. 

(1) In section[3]we present numerical simulations of the nonlinear periodic Maxwell equations, (1.4), for 



initial data obtained from the explicit NLCME soliton. Under this time-evolution there is robust 
spatially localized structure on the scale of the NLCME soliton envelope. The persistence of a 
localized structure and speed of propagation are consistent with that of the NLCME soliton. There 
is, however, a deviation from the NLCME soliton related to third harmonic generation; these are the 
two accessory pulses around the principle wave in Figure [2] (a) . 
(2) On the microscopic scale of the carrier there is nonlinear steepening and shock formation. Therefore, 
the solution does not evolve as a slowly varying envelope of a single frequency carrier wave. The long- 
lived and spatially localized coherent structure which emerges has the character of a slowly varying 
envelope of a train of shocks. We call this an envelope carrier-shock train. Figure [3] illustrates the 
shock-like small spatial scale behavior under slowly varying envelope. 



(3) Numerical solution of the nonlinear Maxwell equation (1.4) is non-trivial due to the cubic nonlin- 
earity. As a hyperbolic system, it is neither genuinely nonlinear nor linearly degenerate 23 29 44 . 
To solve by finite volume methods, as we do, an explicit solution of the Riemann problem must be 
constructed. Details of this are given in Appendix [B"} 

The appropriate entropy condition could, in principle be derived from physical regularization 
mechanisms, which play the role of viscosity in gas dynamics. However, these mechanisms are not 
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Figure 3. On the left is a simulation of the Maxwell equations. On the right is the 
simulation of a truncated asymptotic system. Both simulations were initiated with the same 
initial conditions. There is an indication of shock formation in the left. On the right, we 
see that once sufficiently many harmonics are included, the Gibbs effect appears, confirming 
shock formation. 

well understood. However, such mechanisms and the appropriate notion of weak solution would 
respect thermodynamic principles, which are built into our numerical scheme. 
(4) Using a nonlinear geometric optics expansion 6,13 16,21,30 , systematically keeping all resonances, 
we obtain nonlocal equations governing the interaction of all forward and backward propagating 
modes. Our asymptotic nonlocal system captures the slowly varying envelope of carrier-shock struc- 
tures described above; see below. 



(1.7) 



(1.8a) 



(1.8b) 



(1.9) 



Specifically, we introduce the general wave-form (much more general than a slowly varying enve- 
lope of a nearly monochromatic carrier plane wave), which includes all harmonics 



E(z,t) = e 1/2 ( E+(z-t,ez,et) + E~ (z + t,ez,et) + 0(e)) 



Let 



cf>± = z =p t, et = T, and ez — Z . 

At leading order, the slow evolution of backward and forward components is governed by the coupled 
integro- differential equations: 



d T E+ + d z E+ = dj, (N(<f> A 



^ + 2s,Z,T)) 

1 



s)E~ 

*3 



d T E- - d z E- = -d 4 (N{cfi_ - s)E+{cf>- - 2s, Z, T)). 



-Td d> 



E' 



E- 



Here (•) is an averaging operation in the <j) argument; 



( / ) = Jim 



1 



/(s) ds; 



see also section [4] Equations ( 1.8 1 arise as a constraint on E ± (4>± , Z, T) ensuring that the 0(e) error 
term (1.7) in remains small on large time scales: T = 0(1) or equivalently t = 0(er 1 ). 



Spatial variations in the refractive index, N(z), give rise to a coupling of backward and forward 



waves. Indeed, if N(z) = and one specifies data for the system ( 1.8 1 at t = with non-zero forward 
components (E + ^ 0) and, no backward components (E~ 

5 



0) then, formally, E remains zero for 
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Figure 4. Truncating (1.10) to odd harmonics \p\ < 16, we simulate the initial value 
problem an NLCME soliton in the first harmonic, and the others zero. The above time 
series of the energy associated with each harmonic, e p , shows that most of the energy 
continues to reside in the first harmonic. 



all time, i.e. no backward waves are generated. Continuing with this assumption of N = and 
Eq = 0, if we let V(</>, T) = E+{<j), Z — T,T), with Z„ arbitrary, then V satisfies 

This generalized Burger's equation will gives rise to a finite time singularity. We revisit this obser- 
vation in the discussion, section [6j when considering how singularities might appear when the linear 
coupling between backwards and forwards waves is restored, N 0. 
(5) The nonlocal equations may also be written as an infinite system of coupled mode equations. In the 



case where E is 2ir— periodic in cf>±, the integro-differential equation (1.8) reduces to an infinite 



(1.10a) 



system of coupled mode equations for the Fourier coefficients {Ep(Z, T) :peZ} 
d T E+ + d z E; = i P N 2p E~ + \p T - E+E+E^_ q _ 



(1.10b) 



d T E- - d z E- = i P N 2p E+ + ip^ E-E-E, 



p-q- 



We call this system the extended nonlinear coupled mode equations (xNLCME). xNLCME reduces 
to the classical NLCME if we neglect higher harmonics. 



(6) Simulations of successively higher dimensional mode truncations of (1.10) show improved resolution 
of the carrier shocks under a slowly varying envelope, whose scale is captured by a comparatively low 
order truncation. Indeed, Figure [2] (b) shows that inclusion of the third harmonic in the asymptotic 
system resolves the large scale feature, while inclusion of additional harmonics in Figure [3] (b) shows 
the Gibbs effect, expected for a finite Fourier representation of a discontinuous function. This 
demonstrates that our asymptotic analysis leads to equations capturing the essential features of 
nonlinear Maxwell. However, if we consider how energy, initially only in the first harmonic, is 
redistributed in time, we see in Figure [4] that most of the energy persists in the first harmonic. This 
reflects the partial success of NLCME as a model for periodic nonlinear Maxwell. 



arc 



Relation to previous work: Some of the earliest examinations on optical shocks can be found in Rosen, 39 
and, DeMartini et al. [5] . In these works, the authors applied the method of characteristics to a unidicretional 
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model. Kinsler and Kinsler et al. have continued to examine this problem, and have developed an algorithm 



for detecting the onset of shock formation, 18 19 . Carrier shocks were also examined by Flesch, Moloney, 
& Mlejnek, [9j, for spatially homogeneous Maxwell system with chromatic dispersion, modeled via a time- 
nonlocal Lorentzian polarization response. Ranka, Windeler, & Stentz have found experimental evidence of 
optical shocks, |37| . In their work, a monochromatic pulse with sufficient power steepened and generated a 
broadband optical continuum. 

Coherent structures in nonlinear and periodic media have also been studied by LeVeque, LeVeque & 



Yong, and Ketcheson 17 25 26 in a model for heterogeneous nonlinear elastic media. They considered 
order one solutions in high contrast, rapidly varying, periodic structures. Their simulations yielded localized 
structures on the scale of many periods with oscillations on the scale of the period. For piecewise constant 
(discontinuous) periodic structures, they have a discontinuous carrier shock-like character on the scale of 
the period, though this is due to discontinuities in the medium, the fluxes remain continuous. A two- 
scale (homogenization) expansion yields a nonlinear dispersive equation, with solitary waves, similar to the 
computed solution envelope. In their physical regime, the variations in the properties of the media and the 
nonlinearity are 0(1). In contrast, we consider an asymptotic regime where the constrast of the periodic 
structure and nonlinearity are of the same order, 0(e). Furthermore, the initial condition has two scales 
(envelope and carrier scales), where the carrier wave length is of the same order, indeed in resonance with, 
the periodic structure. These different scalings lead to different asymptotic descriptions. An early example 
of the interactions between nonlinearity and a periodic structure was in atmospheric science, studied by 



Majda et al., 31 . In this work, a model of the interaction of equatorial waves with topography gives rise to 
nonsmooth profiles (in this case, solitary waves with corner singularities). 

Finally, systems of coupled modes have also been examined in prior works, though the work is typically 
limited two just two harmonics, such as a first and second harmonic system or a first and third harmonic sys- 



tem. Such a system was studied by Tasgal, Band, & Malomed 46 , who were able to find stable polychromatic 
solitons in a first and third harmonic system. 

An outline of this paper is as follows. In Section [2j we review how NLCME arises as an approximation of 
nonlinear Maxwell. Results of Maxwell Simulations, showing the coherent structures and shocks, are given 
presented in Section |3]VVe then present our derivation of xNLCME in Section [4j followed by simulations of 
this system in Section[5j We discuss all of these results in Section [6] 

Acknowledgements: The authors would like to thank R.R. Rosales for discussions during the early 
stages of this work on the use nonlinear geometrical optics. We also thank M. Pugh, D. Ketcheson, R.J. 
LeVeque, and C. Sulem for helpful discussions. GS was supported in part by NSF-IGERT grant DGE-02- 
21041, NSF-CMG grant DMS-05-30853, and NSERC. MIW was supported in part by NSF grants DMS-07- 
07850 and DMS-10-08855. MIW would also like to acknowledge the hospitality of the Courant Institute of 
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2. Nonlinear Maxwell and NLCME 

In this section we briefly review how NLCME arises from nonlinear Maxwell with a periodically varying 
index of refraction. We also identify the mathematical inconsistency of NLCME as a description of the 
wave-envelope. 



First, we write the nonlinear Maxwell equation (1.4 1 as 

(2.1) d 2 t (n(z) 2 E + X E 3 ) = d\E 
with index of refraction 

(2.2) n(z) = 1 + eN{z), < e < 1, 
where N(z) is 2ir periodic with mean zero and Fourier series: 

(2.3) N(z) = J2 N P eipZ - 

P ez\{o} 

We shall seek solutions which incorporate (i) slow variations in time and space, due to the weak modulation 
about a constant refractive index; (ii) a scaling of the wave-field which seeks solutions in which the effects 



of dispersion and nonlinearity are in balance: 

(2.4) E e (z,t) = £ 3 £ e (z,t;Z,T), Z = ez, T = et. 

Rewriting ( |2.1[ ) in terms of new variables dependent £ e and independent (z, t, Z, T) variables, we obtain: 

(d?-d*)£ € + e(2d t d T £ e -2d z d z £ £ + 2N(z)£ e + X (£ c f ) + 0{e 3 ) = 0. 
Formally expanding £ 6 as 

£ e {z,t,Z,T) = S (z,t,Z,T) + tS 1 {z,t,Z,T) + ... 
we obtain the following hierarchy for £j(z, t, Z, T), j > 0: 
0(6°) {dl-dl)£^Q 

(Die 1 ) (d 2 t ~ d 2 z ) £ x = -2d t d T £ + 2d z d z £ - 2N(z)£ - x (£ Q f 



(2.5) 



0{eP) (d 2 — d z ) £j = expressions in terms of < I < j — 1 



Solving the O(e ) equation yields: 
(2.6) £ (z,t,Z,T) =£+(Z,T)e l{z - t) + £'(Z,T)e' l( - z+t ^ + c.c. 

Thus, the leading order consists of backward and forward propagating waves, modulated by the slow envelope 
amplitude functions £ ± (Z,T), which are to be determined. 

Substitution of ( |2.6| ) into the 0(e 1 ) equation for £i yields the equation: 

(df - dl) £ X 





2id T £ + 


- 2id z £ + 


- 2N 2 £~ 


-3x( 


[\£+\ 2 - 


K2|£-| 2 ) £+ 




(2.7) 

+ 


2id T £~ 


- 2id z £ + 


- 2N 2 £+ 


-3x( 


>"| 2 H 


-2\£+\ 2 )£- 





+ (£ +) 3 e 3l(2 - t} + (£-) 3 e ~ Mz+t) + c.c. + non-resonant terms 
We have used that N = and 

N{z){£ + e l(z - l) + £- e- l(z+t] ) 

= N_ 2 £ + er l{ - z+t) + N 2 £~e i[z ~ t] + c.c. + non-resonant terms. 



(2.8) 



Each term, explicitly written on the right hand side of (2.7), is resonant with the kernel of (df — . It 
follows that the coefficients of all harmonic plane waves: e ±lql(Z ~ t ^ and e ±l i( z + t ) ^ q £ Z must vanish for £\ 
to be bounded in t. 

The vanishing of the coefficients of e^ z ~*^ and e~ l ( z+t ) yields the nonlinear coupled mode equations 
(NLCME): 



(2.9a) 
(2.9b) 



d T £ + + dz£ + = iN 2 £~ 
d T £~ - d z £~ = iN 2 £^ 



iTl \£+\ 2 + 2\£-\ 2 )£ + 



■iTl\£-\ 2 + 2\£ + \ 2 )£- 



ll 



NLCME also has 



where T = |x an d N 2 = The initial value problem for (2.9 1 is well-posed 

explicit family of gap-soliton solutions; see Appendix [XJ 

However, requiring £ ± to satisfy (2.9) removes only the lowest harmonic resonances. This is the approxi- 



mation invoked in the physics literature; see the survey 14] and references cited therein. 



Note however that the remaining explicitly displayed terms on the right hand side of (2.7) are resonant as 



well and induce linear in time growth. If we choose to remove the resonant terms proportional to e 3l ^ _t ^ and 
e -3i(z+t) ky mc i u ding slow modulations of these plane waves at 0(e°), nonlinearity and parametric forcing 
through N(z) will generate yet other resonant harmonics. 



A leading order solution which does not generate resonant terms at higher order must contain all harmon- 
ics. Thus, NLCME is mathematically inconsistent. In section^ we derive an integro- differential equation, 
which consistently incorporates all resonances. As seen from our numerical and asymptotic studies, this non- 
local nonlinear geometrical optics equation more accurately capture features on both small and large spatial 
scales, e.g. changes in the envelope due to higher harmonic generation, as well as carrier shock formation. 



3. Simulations of nonlinear periodic Maxwell 
In this section we discuss the results of numerical simulations, based on the algorithms of Appendix [Bj of 



the nonlinear and periodic Maxwell equations (2.1 1. 

• In section [3TT| we show that for Cauchy initial data derived from the classical NLCME soliton, there 
evolve spatially localized soliton-like states which persist on long time scales. We discuss aspects of 
the large scale (envelope) structure of such states, which are consistent with the NLCME soliton, as 
well as significant deviations. 



In section 3.2 we show that smoothness breaks down in finite time. In particular, we observe shock 
formation on the fast spatial scale of the carrier wave, while a slowly varying envelope evolves 
smoothly. 



We begin by expressing (2.1 ) as a first order system 



(3.D 3 ,("« 2 v x£ >a, (:*)=». 

We introduce the scaling (E, B, D) T = e 1/2 {E, B, D), and expressing the equations in terms of the variables: 
(D, B) coordinates. Dropping tildes, this is 

where E(D, z) is the unique real solution of 

(3.3) D = n{z) 2 E + e X E 3 



3.1. Soliton-like coherent structures. As is well known [Tj,[3J NLCME has spatially localized gap soliton 
solutions. We use the analytical expression for the gap soliton to generate Cauchy initial data, E(z, 0), d t E(z, 0) 



for (2.1) and numerically simulate the evolution. 

Using (2.6) and the leading order approximation for the magnetic field Bf = ^Bj , NLCME soliton data 



(see (A.l) in Appendix | A[) can be seeded into Maxwell using 

(3.4a) E = £+(ez, etje*** - *) + £~{ez, et)e- 4 (*+*> + c.c., 

(3.4b) B = -£+(ez, etje'** - *) + £+(ez, tt)e 1 ^ + c.c. 



We obtain D via (3.3) and evaluate at t = to get the initial condition. 
For a spatially varying index of refraction, we take 

(3.5) N(z) = |cos(2z), i.e. N 2 = N_ 2 = f, N p = 0, \p\ ^ 2, 

e = 0.0625, and x = 1 (L = §). The results of our simulations appear in a - d of Figures |]and[6j While there 
is attenuation in amplitude and some dispersive spreading of energy, the solution remains spatially localized 
over long time intervals. Not only is there a persistence of the localization (with the periodic medium), but 
also there is good pointwise agreement with the NLCME approximation; see Figure [7| 

Frames e - h of Figures [5] and [6] display the corresponding results in the absence of a periodic structure, 
i.e. N(z) = 0. The derealization, dispersive spreading and attenuation of the wave amplitude is greatly 
enhanced. To understand this heuristically, note that a gap soliton is a localized state whose frequency lies in 
the spectral gap of the linearized PDE at the zero solution. A focusing nonlinearity adds a (self-consistent) 
potential well, creating a (nonlinear) defect mode with frequency lying in this spectral gap. If N(z) = then 
the linearization at the zero state has no spectral gap. Thus, a oscillating with the gap soliton frequency 
would couple to radiation modes and dispersively spread and attenuate. This mechanism is discussed, for 



example, in 45 



Simulations with varying refractive index, (3.5) 
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(c) (d) 
Simulations with constant refractive index, N(z) = 0: 
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Figure 5. Rescaled Maxwell equation, (3.!?), time-evolution for data generated by the 
NLCME soliton with parameters v = .9 and S — .9; sec (A.ll. The solutions are computed 
with 20000 grid points on the domain [-500, 500]. 



Simulations with varying refractive index, (3.5) 
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(c) (d) 
Simulations with constant refractive index, N(z) = 0: 
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FIGURE 6. Solution of rescaled nonlinear periodic Maxwell equation, ( |3.2[ ), for initial data 
generated by the NLCME soliton with parameters v — and S — 7r/2; see (A.ll.The 
solutions are computed with 20000 grid points on the domain [—500, 500]. 
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Figure 7. Comparison of the solution appearing in Figure [5] a - d, with the exact NLCME soliton. 



We note that it is also essential that the data be properly prepared to see a persistence of localization. 
For the initial condition 



(3.6a) 
(3.6b) 



D = 0.5cos(z)sech(e2), 
B = -D, 



we see in Figure [H] substantial spreading. This data mimics the gap soliton's amplitude, slowly varying 
envelope, and carrier wave, but is apparently too far outside the basin of attraction to converge to a localized 
state. Similar results were observed with Gaussian wave packet initial conditions. 



3.2. Envelope carrier-shock trains. Although ther the slowly varying NLCME envelope shape is robust, 
for the nonlinear Maxwell time-evolution, there is evidence of nonlinear steepening and shock formation on 
the short (carrier) microstructure spatial scale. Thus, the nearly monochromatic slowly varying envelope 
approximation of NLCME is violated. 

Figure [9] displays the time-evolution for (a) moving and (b) stationary NLCME - gap soliton data. For 
each initial condition, the nonlinear Maxwell evolution is simulated for different grid spacings. As we increase 
the number of grid points, sharp features are better resolved by the shock capturing algorithm. One can 
also examine the Fourier transform of the output to see that we obtain an algebraically decaying solution in 
wave number, with peaks at the integer wave number values. 

In summary our observations support the emergence of an envelope carrier-shock train; persistence of 
coherent, slowly varying, wave envelope and shock formation on the carrier scale. 
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Figure 8. Solution of rescaled nonlinear periodic Maxwell equation, (3.2) with periodic 



refractive index (13. 5|), for initial data (3.6) . In contrast to the NLCME soliton data, the 



shape of the solution does not persist. The solution is computed with 20000 grid points on 
the domain [-500,500]. 
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Figure 9. Increasing the number of grids points better resolves the shocks in the carrier 



wave. For NLCME soliton data with the indicated v and S (see (A.l)), with index of 
refraction N(z) given by (3.5)). 
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4. Resonant nonlinear geometrical optics and nonlinear spatially inhomogeneous 

Maxwell equations 

In this section we derive a system of equations, which incorporates all wave-resonances and which our 
numerical simulations show, captures the key features of the nonlinear Maxwell time-evolution, in par- 
ticular, the presence of robust envelope carrier-shock train solutions. We derive this system, for general 
non-homogeneous media, using a nonlinear geometrical optics expansion; see, for example, 13 14 ,30 . The 
equations obtained are the general integro-differcntial equations (1.8). In the case of a periodic medium, 
they reduce to an infinite set of local equations, which we call the extended nonlinear coupled mode equations 
(xNLCME). If, in xNLCME, we neglect all but principle resonances, xNLCME reduces to NLCME. 

As we shall see, in our numerical simulations of increasing high dimensional truncations of xNLCME (sec- 
tion [5]), this theory appears to accommodate the observed carrier shocks and large scale coherent structures. 

4.1. Nonlinear geometric optics expansion. In contrast to the ansatz of Section [2j we assume the more 
general form 

(4.1) u(z, t) = u (0) (z, t, Z, T) + eu (1) (z, t, Z, T) + e 2 u (2) (z, t,Z,T) + 

where u = (E, B) T and Z = ez, T = et. Inserting (|4.1[) into (|3.2|, (|3.3[), the first order system 



dt 



n(z) 2 E + exE 3 
B 



-B 
-E 



we expand to get, 



(d t + B^d z ) u<°> + e \{d t + B^d z ) 



with matrices 
(4.2) 

At O(e ), 
(4.3) 



B (0) = 



-1 



+A ( - 1 \z, u)S t u(°) 



(d T + B^d z 



,(o) 



2N(z) + 3x^ 2 




0. 



Solving this as the generalized Eigenvalue problem, 



XI r = 



the solutions are: 
(4.4) 

The corresponding left eigenvectors are 



A±=±l, r± = 



(i 



(4.5) 

With this normalization, liA^'Tj = Sij 
(4.6a) 
(4.6b) 
(4.6c) 



1± = |(1 Tl) 
The leading order fields are then 

(°> = E+ (</>+, Z,T)v+ +E- (</>-, Z,T)t- 



E<-V = E+ (0+, Z, T) + E~ (<f>- , Z, T) , 
4>± = z T t. 



This expression is much more general than (2.6) used in the derivation of NLCME. 
At 0(e), the equation is 



(4.7) 

If we assume 
(4.8) 



dt 



B^d z j u« =-(dr + B^d z j u<°> - A^(z, u^)d t u^ 
u*- 1 ' (z,t) = m + (z, t)r + + m~(z, t)r_, 



and substitute into (4.7), then left multiply by 1 + and then by 1 , we get the two equations 

- (d t m+ +d z m+) = d T E+ +d z E+ + l + A^{u^)x 
(4.9) y-u^ x + 



- (d t m- - d z mT) = &tE~ - d z E- + l_A (1) (u (0) )x 
(4.10) (-d lP+ E+r+ + d^_E~vJ) . 

The last term is the same in both equations, 

l±4 (1) (u (0) ) (~d <p+ E+r+ + fy_£Tr_) 

= i (2N(z) + 3 X i? (0)2 ) (-d^ + E+ + d^E-) 
= N(z) {-d^ + E+ + d^_E-) 



(4.11) 



|X (E+ + E-y {-d < p + E+ + d+_ET 



Integration of (4.9) along the characteristic dtz + = 1 from t = to t = L, yields 

-(m+(z+(L),L)-m+(z + (0),0)) = 

/ d T E+(Z, T, z+(0)) + d z E + (Z, T, z+(0))ds 
Jo 



(4.12) 



JV(z + (s))c> 0+ £ + (Z,T, z+(0))ds 



+ / N(z + (s))d <t> _E-{Z,T,z + (s) + s)ds 



\ x (E+(Z, T, z+(0)) - E-(Z, T, z+(a) + a)) d 



xd^ + E+(Z,T, z+(0))] ds 



X (E+(Z, T, z+(0)) - E-(Z, T, z + {s) + a))' 



xd 4> _E-{Z 1 T 1 z + {s) + s)] ds. 
Similarly, integration of ( 4. 10[ ) along the characteristic <9tz_ = — 1, yields 

-(m-(z_(L),L)-™-(z_(0),0)) = 

/ d T E-(Z, T, z+(0)) - d z E- (Z, T, z+(0))ds 
Jo 



(4.13) 



iV(z_(s))d 0+ £ + (Z, T, z-(s) - s)ds 



+ / N(z_(s))d 4> _E- (Z,T, z_(Q))ds 



L r 



X (E+(Z,T,z-(s) -a)- E~(Z,T, z_(0)))' 
d^ + E + (Z,T, z-(a)-a)] ds 



X (E+(Z,T,z-(a) -a)- E~(Z,T, z_(0))) J 



:d c j > _E-(Z,T, z_(0))] ds. 
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Necessary conditions for m± to grow sublinearly in t as t — > oo are the solvability conditions: 
d T E+ (Z,T,z+(0)) + d z E+(Z, T, z+ (0) ) = 
1 



(4.14a) 



(4.14b) 



(4.14c) 



lim 

L— >oo Li 

lim — 



N{z + {s))d z+[Q) E-{Z, T, z+(s) + s)ds 
|X (E+(Z, T, z+(0)) + E-(Z, T, z+(s) + a)f 



xd z+{o) E+(Z,T,z + (0))] ds, 



8 T E- (Z, T, z_ (0)) - d z E- (Z, T, z_ (0)) = 



lim — 

L—too L 



lim — 

L oo L 



N(z-(s))d z _ i0) E+(Z, T, z_(a) - s)ds 



X (E+(Z, T, z_ (s) -s)-E~ (Z, T, z_ (0)))' 



d 



z_(0) 



E~(Z,T, 2_(0))] ds. 
Given (z,t), z+(0) = z — t = and z—(0) = z + t = 

(4-15) (/) 



lim — 



_. Defining 

L 

f(s)ds 



the equations may be compactly expressed as: 

d T E+ + d z E+ = - (N((f>+ + 8)8^- 
[[E+f + 2E+ (E~) + ((E- 



(4.16a) 
(4.16b) 



fx 



d^E+, 



d T E- - d z E- = (N(0_ - s)d 4> E + {(j)- - 2s)) s 
-l X ((E-) 2 + 2E-(E+) + ((E+) 2 ))d,E- 

are d 



It is important to recognize that the arguments of the fields in ( |4.16a 
(4.16b), they are 



= z — t, Z, and T, while in 



z + t, Z, and T. As in our derivation of NLCME in Section [2j T = |%. With this 
notation, (4.16) can be rewritten, after an integration by parts, in conservation law form, 

d T E+ + d z E+ = <9 <ni(0 + + s)E-{<t>+ + 2s) g 

(4.17a) 



+ Yd, i (E+) 3 + (E+f (E-) + E+ ((E-) 2 
d T E- - d z E- = -8+ (m(0_ - s)E+{4>_ - 2s)) 



(4.17b) 



E~ 



(E+) (E- 



3 v ; 



Equations ( |4.17[ ) corresponds to the integro-differential equations of the introduction, if we omit the (E^) 
terms. Since (E =t ) is time-invariant (see section 4.3) by choosing initial conditions for which (E^) (T = 0) = 
0, these terms can be dropped from (4.17). Finally, note that (4.16) are applicable to a general heterogeneous 
dielectric material with the appropriate scalings. 

4.2. Periodic Media and xNLCME. We now specialize to the periodic case. Assume now that N(z + 
2ir) = N(z). Then (4.17) is invariant under the discrete translation: <j> H> <fi + 2ir, i.e. 

E+((f>, Z, T) ^ E+{<f> + 2ir,Z, T) 



(4.18a) 
(4.18b) 



E~(4>,Z,T) i y E~ 



2tt, Z, T) . 



Thus, under the assumption of existence and uniqueness of solutions to (4.17), if the initial data are 2n in 
the <p argument, then the solutions remain 2-7T periodic in <f>. In the periodic setting, the averaging operator, 
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(4.15 ), simplifies to 



(/)= 2?r y o Mda. 
We now expand N(z) and E ± in Fourier series, 

(4.19) N(z) = J2 N P j PZ , 

(4.20) E ± ( ( j ) ,Z,T) = J2E±(Z,T)e ±i ^, 

p 

where N p — N_ p and = E_ since N and E^ are real valued. In this case, the system of Fourier 
coefficients {E p (Z,T) : p G Z} satisfy the infinite system of extended nonlinear coupled mode equations 
(xNLCME): 



8 T E+ + d z E+ = i P N 2p E- + ip- 



(4.21a) 



Y,E+E+E+_ q _ r 



drB- - = i P N 2p E+ + ip 



+3E» KK- q + 3 ( E K~ I 2 ) £ p 

r 



(4.21b) 



p—q—r 



+3£ + £ + 3 ^ \E+\ 2 )E p 



4.3. Conservation Laws and Hamiltonian Structure. Equation (4.17), and alternatively (4.21), have 
two conservation laws: 



Proposition 4.1. Assume that E ± is a sufficiently smooth and sufficiently Z— decaying solution of (4.17) 
and that {E p (Z,T)} p ez is the corresponding solution of xNLCME. Then, 



(4.22a) 
(4.22b) 
(4.22c) 
(4.22d) 



d 

dT 
d 

dT 



= 



J (E + (;T)) dZ= A J E+(;T)dZ 

J (E+(; T)) dZ=^ J B "(.,T) dZ = 
J ((E+) 2 (;T)) + ((E-f(;T)) dZ 
J \E+(;T)\ 2 + \E-{;T)\ 2 dZ = 0. 



_d_ 
dT 
d 

dT 



Proof. Setting p = in (4.21), 



d T E+ + d z E+ = 0, 
drEo - d z E^ = 0. 



Integrating in Z establishes the first two conservation laws in terms of the Fourier modes. Integrating (4.20) 
in <f> over [0,2tt) relates (E ± ) to Ejf. 
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Multiplying (4.21a) by E£, summing over p, and adding its complex conjugate, 

E E q E r E -p E p-q-r 



Ip 



J2 9t \E+\ 2 + d z \E+\ 2 = WN2 P E~E+ + 3 E 

p p 

+ 3E E E q E P -q E - P 

q 

+ 3 (ei^i 2 ) \k\ 

The quartic terms will all vanish. Consider the first quartic term, and note that 



p—q—r 



p,q,r 



(4.23) 



E ^E+E+E+E, 

k 1 + k 2 +k 3 +k i =a 

E 

k 1 + k 2 +k 3 +k 4 =0 

E 

k 1 + k 2 +k 3 +k i =0 

E ^KKKK 

k 1 + k 2 +k 3 +k i =0 



k 2 Et Ef Et Ef 

* k\ K2 K3 K4 



ksE+E+E+El 



Hence, 



(4.24) 



-q—r 



J2 pE+E+E+ p E+ 

7 E (k 1 + k 2 + k 3 + k^E+E+E+El = 



fci+fc 2 +*:3 + fe4=0 

The second quartic term vanishes using a similar analysis. The last quartic term, 



(4.25) 



EWEK 



will vanish because the p and —p terms will cancel one another. Similar analysis holds for (4.21b I, leaving 
us with the two equations 



(4.26) 


E^KI 2 


+ dz\E+\ 


(4.27) 


E^-f 


-d z \E~\ 


Summing 


these two, and integrating 


in Z gives 



p ' 



he L 2 conservation law. □ 
To simplify our analysis we assume Eq are initially zero from here on. The equations reduce 



to 



(4.28a) 



(4.28b) 



d T E+ + d z E+ = 9 (N(<t>+ + s)E-(<j>+ + 2s) g 

+ Td^(E+) 3 + E+((E-) 2 
d T E- - d z E- = -d+ (N(<t>- - s)E + (<p_ - 2s)) s 

-rd,[((E+) 2 ) E - + l(E-f 

is 



and 
(4.29a) 

(4.29b) 



d T E+ + d z E+ = ipN 2p E p 



■VP-. 



[E 



77*+ rp+ 771+ 
q r p—q — V 



d T E~ - d z E- 



i P N 2p E; + ip^ [J2 E~E-E~_ q 



£,1 



—q—r 
2\ 



E~ 



These are equations (1.8| and (1.10) from the introduction. Truncating (4.29) to just mode recovers 
the NLCME, subject to the identification of £ ± with Ef. 

Another time-invariant functional is a consequence of the Hamiltonian structure given in the following 
result, which is straightforward to verify: 



Proposition 4.2. The system (4.29) is a Hamiltonian system 
(4.30) 



d T E^ = -ip^-, d T E =-ip^= 

se; SE p 



where with time-invariant Hamiltonian 



H[E ± ,E ± ] 



n(-,T) 



and Hamiltonian density 

1 

H(Z,T) = i J2 ~ (KdzEp! ~ E v^E pi ) - E N *n E l E n 



(4.31) 



2 ^ pi 

Pi=i ^ 1 

r 1 1 

3 24 



E 

>1+P2 +P3+P4=0 



Pl=l 

E P 1 E i E P 3 E P4 + ^Pl S P2 ^Pa ^P4 



,11 

22 



\ pi / \ pi / 



c.c. 



5. Simulations of the Truncated xNLCME 

In this section we simulate truncations of the infinite dimensional xNLCME system, performed pseudo- 
spectrally with fourth order Runge-Kutta time stepping. These simulations suggest that 

• xNLCME has its own localized soliton-like structures which better capture the dynamics of the 
nonlinear periodic Maxwell equation for our class of initial conditions than NLCME and 

• xNLCME has singular solutions, {E p (Z, T)} with a cascade of energy to higher wave numbers, p. 
The physical electric field 

E(z,t) « (h (E + (z -t,,ez,et) + E~ (z + 1, ,ez,et) ) 

= ek E (E+(Z,T)e l ^ z - T ^ e + E;(Z,T)e-^ z + T ^ e 
pez\o 
+ c.c. 

develops a carrier-shock train structure. 

As we saw in Section 3T particularly Figure |6j though the NLCME soliton data appeared robust, there 
was some escape of energy. This can be accounted for in the xNLCME through the inclusion of additional 
modes. 

Starting with the same initial conditions, we simulate the NLCME soliton of E^ with soliton parameters 
v = and <5 = £ , and material parameters 

2 



1, N : 



±2 



N 







19 



-500 -400 -300 -200 -100 100 200 300 400 500 -500 ^100 -300 -200 -100 100 200 300 400 500 




-500 -400 -300 -200 -100 100 200 300 400 500 



-500 -400 -300 -200 -100 100 200 300 400 500 



FIGURE 10. Evolution of an NLCME soliton in the xNLCME, resolving odd modes \p\ < 4. 
Computed with 4096 grid points in the Z coordinate. Compare with Figure [6] 



in (4.291 resolving only a finite number of harmonics. The primitive electric field is reconstructed from these 



simulations as 



(5.1) 



E = 



P=-Pm, 



E+{Z,T)M Z - T V* + E-(Z,T)e-^ z+T ^ t + 



c.c. 



E is plotted in Figures [lOj and [TT] which resolve odd modes up to 3 and 15, respectively. Comparing with 
Figure |6j we infer that the two smaller pulses s ymm etrically expelled from the main wave were transferred 
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This addresses the macroscopic discrepancy between 



into E ±3 , since these clearly appear in Figure 
NLCME and Maxwell. 

Including the additional modes also suggests shock formation by re-examining Figure [9] The sharp er, 
shock like features, can only be resolved by the inclusion of the the higher harmonics. The contrast between 
different truncations is shown in Figure [l2j Indeed, we see the Gibbs phenomenon that would be expected 
from taking a truncated Fourier representation of a discontinuous function. 

Despite this, NLCME still gets certain leading order effects correct, such as the main structure in the 
Maxwell simulations. The robustness of NLCME can also be seen by exploring how energy is partitioned 
amongst the harmonics. Let 



(5.2) 



\K\ 



E^ 



e: 



dZ, p = 1, 3, 



This is the energy associated with mode p. Their sum is conserved. Plotting this for the above simulations 



in Figure 13 we see that most of the energy remains in mode one, some migrates into mode three, and less 
in the subsequent modes. 
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FIGURE 1 1 . Evolution of an NLCME solition in the xNLCME, resolving odd modes \p\ < 16. 
Computed with 16384 grid points in the Z coordinate. Compare with Figure [6| 



6. Summary and discussion 

We first numerically simulated the one-dimensional nonlinear Maxwell equations in the regime of weak 
nonlinearity, low contrast periodic structure (weak dispersion) with wave-packet data satisfying a Bragg 
resonance condition, i.e. carrier wavelength equal to twice the medium periodicity. We observe strong ev- 
idence of the emergence of a coherent structure evolving as slowly varying envelope with a carrier-shock 
train. This violates the nearly-monochromatic assumption underlying the classical nonlinear coupled mode 
equations. We propose our nonlocal integro-differential equations governing coupled forward and backward 
waves, derived via a nonlinear geometrical optics expansion, as the physically correct, mathematically con- 
sistent description of waves governed by nonlinear Maxwell in a periodic structure with negligible chromatic 
(nonlocal in time) dispersion. These equations are equivalent to an infinite dimensional system of couple first 
order PDEs, the extended coupled mode system (xNLCME). The electric field, E, obtained from numerical 
solution of successively higher truncations of xNLCME converges toward the envelope carrier-shock trains 
observed in direct simulations of the nonlinear Maxwell equations. 

Finally we mention that our methods could be applied to study the long time evolution of wave-packet 



type initial conditions for the problem of quadratically nonlinear elastic media, consider in 17 25 26 We 



obtain nonlocal equations of resonant nonlinear geometrical optics (or equivalently an infinite family of 



nonlinear coupled mode equations), governing interacting forward and backward propagating waves 42 . A 
difference between the quadratic and cubic case is that the smallest truncated system that retains nonlinear 
interactions contains four modes, p = ±1, ±2. Nonlinear effects occur through second harmonic generation, 
a process well-known in nonlinear optics. 

Open problems and conjectures: As our simulations show, there is agreement between finite mode 
truncations of the integro-differential equations and the primitive Maxwell system. Assessing, and proving 
the time of validity of this approximation is one upon problem. 



(a) Resolves odd harmonics |p| < 2 



(b) Resolves odd harmonics |p| < 4 





(c) Resolves odd harmonics \p\ < 8 



(d) Resolves odd harmonics \p\ < 16 



Figure 12. Comparison of the features that develop on the scale of the medium in different 
truncations of the equations. Including additional harmonics better captures the shocks seen 
in Figure [9] 



Following up on assessing the time of validity, there is also the question of the time of existence and 
the well-posedness of the equations. We expect that solutions of xNLCME for initial data having a finite 
number of nonzero mode amplitudes, e.g. NLCME gap soliton data, will give rise to solutions of xNLCME 
that develop singularities in finite time. The nature of this blowup is expected to occur via a cascade 
to high mode amplitudes (higher index, p), corresponding to modes necessary to resolve the carrier shock 
structure in the small scale. As we mentioned in the discussion, there is clearly singularity formation when 
the heterogeneity is turned off (N — 0), and either E + or E~ is initially zero. It is an open problem as to 
whether this particular mechanism for singularity formation will persist when coupling is restored. 

As pointed out in the introduction, the success in modeling experiments with NLCME suggests that, 
although there is such a (weakly turbulent) cascade, it is only a small part of the optical power that is 
transferred to high wavenumbers and that this energy contributes mainly to resolving the small-scale shocks. 
To explore this, one needs to simulate the xNLCME equations with many more harmonics. Plotting the 
Fourier transform (in the Z coordinate) of the simulations in Section [5] in figure 14 we see that the spectral 
support grows as we increase the number of resolved envelopes (the E~'s). A related question is whether or 
not the primitive Maxwell system, the xNLCME system, or one of its truncations possess genuine coherent 
structures. In [46], the authors found such solutions for a first and third harmonic system. This shall be 
further explored in the forthcoming publication, 



3G 



Finally, our computations in Section [3] invoked of a gas-dynamics entropy condition. Such a condition is 
necessary to use finite volume methods. Although thermodynamically consistent, we do not know whether 
this is the correct regularizing mechanism of electrodynamics. 
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(a) Resolves odd harmonics \p\ < 4 



20 40 60 80 100 120 140 160 180 200 



(b) Resolves odd harmonics |p| < 8 




(c) Resolves odd harmonics |p| < 16 



Figure 13. Energy distribution, (5.2), for truncated xNLCME simulations with different 



numbers of harmonics. In all cases, the energy initially residing in mode one tends to stay 
there. 
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Figure 14. Fourier transforms of the solutions to truncations of the xNLCME equations. 
Increasing the number of envelopes expands the support in Fourier space. 



[14] J. Hunter, A. Majda, and R. Rosales. Resonantly interacting weakly nonlinear hyperbolic waves ii: several space variable. 

Studies in Applied Mathematics, 75:187-226, 1986. 
[15] J. Joly, G. Metivier, and J. Rauch. Resonant one dimensional nonlinear geometric optics. Journal of Functional Analysis, 

Jan 1993. 

[16] J. Joly, G. Metivier, and J. Rauch. Global solvability of the anharmonic oscillator model from nonlinear optics. SIAM 

Journal on Mathematical Analysis, Jan 1996. 
[17] D. Ketcheson. High Order Strong Stability Preserving Time Integrators and Numerical Wave Propagation Methods for 

Hyperbolic PDEs. PhD thesis, University of Washington, 2009. 
[18] P. Kinsler. Limits of the unidirectional pulse propagation approximation. Journal of the Optical Society of America B: 

Optical Physics, 24(9):2363-2368, Jan 2007. 
[19] P. Kinsler, S. Radnor, J. Tyrrell, and G. New. Optical carrier wave shocking: Detection and dispersion. Physical Review 

E - Statistical, Nonlinear, and Soft Matter Physics, 75(6), Jan 2007. 
[20] S. Klainerman and A. Majda. Formation of singularities for wave equations including the nonlinear vibrating string. Comm. 

Pure Appl. Math, 33(3):241-263, 1980. 
[21] D. Lannes. Dispersive effects for nonlinear geometrical optics with rectification. Asymptotic Anal, 18(1-2):111— 146, Jan 

1998. 

[22] P. Lax. Development of singularities of solutions of nonlinear hyperbolic partial differential equations. Journal of Mathe- 
matical Physics, 5:611, 1964. 
[23] P. Lax. Hyperbolic Partial Differential Equations. American Mathematical Society, 2006. 
[24] R. LeVeque. CLAWPACK. http://www.amath.washington.edu/ claw/. 

[25] R. LeVeque. Finite volume methods for hyperbolic problems. Cambridge university press, 2002. 

[26] R. LeVeque and D. Yong. Solitary waves in layered nonlinear media. SIAM Journal on Applied Mathematics, pages 1539- 
1560, 2003. 

[27] T. Liu. The Riemann problem for general 2x2 conservation laws. Transactions of the American Mathematical Society, 
pages 89-112, 1974. 

24 



[28; 

[29 
[30 

[31 

[32 

[33 

[34 

[35; 

[36 

pr 

[38; 

[39 

[4o; 

[41 

[42; 

[43. 

[44 
[45; 
[46 

[47; 

[48; 

[19 
[50 



T. Liu. The Riemann problem for general systems of conservation laws. J. Differential Equations, 18(l):218-234, 1975. 
T. Liu. Hyperbolic and viscous conservation laws. Society for Industrial Mathematics, 2000. 

A. Majda and R. Rosales. Resonantly interacting weakly nonlinear hyperbolic waves i: a single space variable. Studies in 
Applied Mathematics, 71:149-179, 1984. 

A. J. Majda, R. Rosales, E. G. Tabak, and C. Turner. Interaction of large-scale equatorial waves and dispersion of kelvin 
waves through topographic resonances. Journal of the Atmospheric Sciences, 56:24, 1999. 

S. Miiller and A. Voss. The riemann problem for the euler equations with nonconvex and nonsmooth equation of state: 
construction of wave curves. SIAM Journal on Scientific Computing, 28(2):651-681, 2006. 

S. Miiller and A. Voss. The Riemann problem for the Euler equations with nonconvex and nonsmooth equation of state: 
construction of wave curves. SIAM Journal on Scientific Computing, 28:651, 2006. 

D. Pelinovsky and G. Schneider. Justification of the coupled-mode approximation for a nonlinear elliptic problem with a 
periodic potential. Applicable Analysis, 86(8):1017-1036, 2007. 

D. Pelinovsky and G. Schneider. Moving gap solitons in periodic potentials. Mathematical Methods in the Applied Sciences, 
31(14):1739-1760, 2008. 

D. Pelinovsky, G. Simpson, and M. Weinstein. Solitons in the nonlinear maxwell equations. In preparation. 

J. Ranka, R. Windeler, and A. Stentz. Optical properties of high-delta air-silica microstructure optical fibers. Optics letters, 

25(ll):796-798, Jan 2000. 

M. Reed and B. Simon. Methods of Modern Mathematical Physics IV. Analysis of Operators. Academic Press, 1978. 

G. Rosen. Electromagnetic shocks and the self-annihilation of intense linearly polarized radiation in an ideal dielectric 

material. Physical Review, 139(2A):A539-A543, Jan 1965. 

G. Schneider. Nonlinear coupled mode dynamics in hyperbolic and parabolic periodically structured spatially extended 
systems. Asymptotic Analysis, 28(2):163-180, 2001. 

G. Schneider and H. Uecker. Existence and stability of modulating pulse solutions in Maxwell's equations describing 
nonlinear optics. Zeitschrift fur Angewandte Mathematik und Physik (ZAMP), 54(4):677-712, 2003. 

G. Simpson and M. Weinstein. Nonlinear geometrical optics and quadratically nonlinear periodic media. Unpublished 
Notes, 2010. 

D. Sjoberg. On uniqueness and continuity for the quasi-linear, bianisotropic Maxwell equations, using an entropy condition. 

Progress In Electromagnetics Research, pages 317-339, 2007. 

J. Smoller. Shock waves and reaction- diffusion equations. Springer, 1994. 

A. Soffer and M. Weinstein. Time dependent resonance theory. Geometric and Functional Analysis (GAFA), 1998. 

R. Tasgal, Y. Band, and B. Malomed. Gap solitons in a medium with third-harmonic generation. Physical Review E - 
Statistical, Nonlinear, and Soft Matter Physics, 72(1):1-10, Jan 2005. 

B. Wendroff. The Riemann problem for materials with nonconvex equations of state. I. Isentropic flow. J. Math. Anal. 
Appl, 38:454-466, 1972. 

B. Wendroff. The Riemann problem for materials with nonconvex equations of state. II. General flow. J. Math. Anal. Appl, 
38:640-658, 1972. 

P. Wesseling and D. Van Der Heul. Uniformly effective numerical methods for hyperbolic systems. Computing, 66(3) :249— 
267, 2001. 

P. Wesseling, D. van der Heul, and C. Vuik. Unified methods for computing compressible and incompressible flows. In 
Proceedings of ECCOMAS, pages 11-14, 2000. 



Appendix A. The NLCME Soliton 



Using the notation of 11 , the NLCME soliton solution of (|2.9| is given by 
(A.la) 



(A.lb) 
(A.lc) 

(A.ld) 
(A.le) 
(A.lf) 



£ (Z, T) = sae 



£-(Z,T) = -ae ir> 




^-sin(5e is,T sech((9- isS/2) 



Asin<5e is<T sech(0 4-is(5/2) 



= -yN 2 smS(Z -vT), a = -yN 2 cos S(vZ - T) 



£ 28 + e -isS\ 2"/(3-^) 
e 20 + e isS 



7 = l/\/l-v 2 , A = 



1-v 
l + v 



s = sign(A^r), a 



2(1 -v 2 ) 



1/4 



K = kNo 



Z-v 2 

We assume that N 2 E M. There are two free parameters, |u| < 1 and 5 E 
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Appendix B. Simulating the Nonlinear Maxwell Equations 







In vector notation, the rescaled Maxwell system, (3.2), and constitutive law, (3.3), are expressed as 

(B.l) *(fl) +fl '(-2^A*) 

d t v + d,f(y, z) = 0. 

To simulate this system of conservation laws, we employ a shock capturing finite volume scheme with the 

Furthermore, we employ the /-wave method to accommodate the spatially 
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varying flux function, [2 

To use finite volume methods we must provide the algorithm with a, possibly approximate, solution of the 
Riemman problem. This introduces a subtlety as our system has a non-convex flux function. Non-convex 
fluxes lead to interesting waves, including rightward (or leftward) traveling rarefaction and Shockwaves 
that are "glued" together. Such waves, sometimes called compound or composite waves, were discussed 
in 



27 28 47 48 and more recently in 33 49 50 . Examples are also give in the texts 25 44 



B.l. Finite Volume Methods for Maxwell. In finite volume numerical methods, at each time step, we 
must solve a Riemann problem between adjacent grid cells: 



(B.2) 



v t + f(v;z,-)a = 
v t + f(v; Zj+i) z = 



for 



< z < z 



.7 + 1/2 > 



for z 



3+1/2 < Z < Zj+3/2, 



v(*,t = t B ) = 



for z 



J- 1/2 



< z < z 



.7 + 1/2 > 



"j+l 



for Zj+l/2 < Z < Zj+3/2- 



z j+i/2 is the interface between the cell centered at Zj and the cell centered at Zj +1 / 2 - The fluxes are assumed 
to be constant in z within each computational cell. We aim to provide an exact solution to the Riemann 
problem, in contrast to an approximate solutions such as the Roe average. 
In the next few sections, we adopt the notation 

v f + f;(v) z = forz<0 

Vt + f r (v) z = for z > 
(B.3) r 

v(z,0) = r 



for z < 
for z > 



where f/(v) = f (v; Zj), f r (v) = f (v; Zj+i), and we take Zj +1 / 2 = 0. 

Given any point v in (D, B) phase space, we construct two, entropy condition (specified below) satisfying, 
manifolds Wi(v) and W2(v). These define the locus of points that can be joined to v by a left-going wave 
in the former case and a right-going wave in the latter case. We parameterize them in the D component. 
Given a state vo, Wj(D;vo) is the parametric curve such that: 

D 

Were the medium homogeneous, solving the Riemann corresponds to finding the state v* that is the 
unique point in Wi(vj) f] W^vv). In terms of the parametric curves, this point solves the equation: 

K W 1 (D i< ;y l ) / 

As the medium is not homogeneous, we match the flux at the interface. We seek and such that: 

(B.5a) W 1 (Df;y l )=Bf 

(B.5b) f i (vf)=f r (v*) 

(B.5c) W 2 (D r ;v*) = B r 

is the the entropy satisfying state immediately to the left of the interface and v* is the entropy satisfying 
state immediately to the right of the interface. 
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(B.4) 



W 2 D 



for D e R } = VV,-(vo) for j = 1, 2 



— B r 



For this problem, the flux matching condition is: 



(B.6a) 
(B.6b) 

Defining transfer function, T, that, 
(B.7) 



E{D* l ;z l ) = E{D*;z r ) 
B\ = B* 

given zi, z r and a left state D^, the flux matched displacement is: 
T{D* l ;z h z r )=D* r 



With this function, (B.5) becomes 
(B.8) 



W 2 D 



B r 



(T{D*- Zl -z r f 

Once we have this value, we recover E* and B* allowing us to compute the fluxes. 
Subject to the specification of the phase space functions Wj are specified, this is a one-dimensional root 
finding problem. 



D* is the unknown 



B.2. Non-Convex Fluxes and the Entropy Condition. It remains to specify the manifolds Wj. This 
requires an additional, non-trivial, assumption on an entropy condition. While such a condition is readily 
apparent in gas dynamics and elasticity, the appropriate condition for Maxwell is non-obvious. 

In this work, we employ a diffusive entropy condition, akin to that found in gas dynamics. This was 
suggested by Sjdberg 43 , as part of an entropy-flux pair involving the Poynting vector. This is also a 
physically consistent, as many dielectrics absorb the higher harmonics that would appear as the wave began 
to shock. 

In constructing the entropy satisfying Wj functions, we closely follow 27 28 32 48 and particularly the 
p-system example in 47 . Graphically, the Wj functions can be constructed by tracing an appropriate convex 
hull of E(D;z). Shock waves occur when points are joined by chords, rarefaction waves when points are 
joined along E(D; z), and composite waves when convex curve is a combination. 

Throughout this section we suppress the z argument, and E'(D) = doE(D,z). 

Since the flux function is no longer uniformly convex, the Lax entropy condition may not be appropriate. 
Instead, Liu entropy condition 27 32 may apply. Recall, that if 

■E(D) 



(B.9) 
then: 



C7(v ,v) 



E(D ) 



B — Bn 



• A shock joining vo to vi satisfies the Lax entropy condition if the system is convex and either 
(B.lOa) Ai(vi) < ct < Ai(v ), a < 
(B.lOb) A 2 (vi) < a < A 2 (v ), a > 

where Ai < < A2 are the eigenvalues of 

-1' 
-E'{D) 

• A shock joining vo to Vi satisfies the Liu entropy condition if 
(B.1I) o-(v ,vi) < cr(v ,v) 

for all points v between the two points along the shock curve in phase space. 

B.2.1. Left Traveling Waves. Given the state Uo = (Dq,Bq) t , we construct Wi(D;uq). Since we have an 
inflection point in E(D) at D = 0, we dissect all the possible configurations of Dq, D and the inflection 
point. Let D\ be the value of D at which the line tangent to (Di, E(Di)) intercepts (Dq,E(Dq)). 

First, suppose that D < Dq < 0, as in Figure |15(a)| In this region, there is no difficulty applying the Lax 
entropy condition (B.IO); there is no shock as 



Ai(-D) = -y/E'(D) < Ai(D„) = -y/E'(D ) < 



Consequently, 



B = B n 



D 
27 



y/E'(s)ds 




Figure 15. Entropy satisfying leftward traveling waves when Dq < 0. 



Still assuming that Dq < 0, if Do < D < D\ the Lax condition continues to apply. Dq and D will satisfy 
(B.10) and we can join the two with a shock, as in Figure |l5(b)| 

B = B + y/\E(D) - E \\D - D \ 

Once D > D\, the solution changes. It is no longer appropriate to apply the Lax condition as we lose 
convexity here. Applying the Liu condition, (B.lll, we see that there is no longer a shock. Indeed, we can 
compute that were there shock solutions, 



o-(vo,v) = 
0"(vo,vi) = 



—E - 


\-E 1 


E — Eo 


B — 


Bo - V 


D-D 


-E x 


+ E Q 


1 Ei — E 


B l - 


-Bo ' V 


Di-D c 



Examining Figure 15(c) 



Ei — E E - E 



implying 



Di - D D - D 
ct(v , vi) < o-(v , v) 



which violates (B.lll. 

As a shock fails to connect the two states, we resort to joining the states by a compound wave. The 
solution is a shock from Do to D\ which continues into a rarefaction wave from D\ to D. Thus, 

cD 

■ "is 



B = B + ^\E(D 1 )-E \\D 1 -Do\+ f ^W^)ds 



At Do = the system is convex yielding leftward traveling rarefaction waves for all values of D: 

B = B + / y/E'(s)ds 

JDo 

For Do > 0, There are again several cases. For D < Di, we have the compound wave again, as in Figure 



16(a) 



B = B - y /\E(D 1 )-E \\D 1 -D \ + / ^{sjds 
As D increases in value and D\ < D < Do, we have a shock solution 

B = B - ^\E(D)-E \\D~Do\ 
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Figure 16. Entropy satisfying leftward traveling waves when Dq > 0. 



See Figure [l6(b)1 Lastly, for D > D , we get a leftward traveling rarefaction: 

f D 

B = B + / y/E'(s)ds 
Jd 

B.2.2. Right Traveling Waves. For right traveling waves, the structure is similar. Given our point Do, let 
(D 2 ,D 2 ) be the point intercepted by the line tangent to (D ,D(D )). 

Again, we first treat Do < 0. The different cases are diagrammed in Figure 17 For D < Do < 0, there is 
a shock, 

B = Do + v / |D(D)-D ||D-Do|. 
For Do < D < 0, this changes to a rarefaction wave, 

i-D 



B = B - [ V / E'(s)ds. 



Crossing the inflection point, < D < D 2 , it becomes a compound wave which rarefacts to the point D* 
followed by a shock, 



/■£>. 

B = B a - ^/W{s)ds - v / |D(D)-D,||D-D.| 

J Dr, 



D* is the point (D*, E±) on the curve whose tangent intercepts (D, E). Past D 2 , the compound wave reduces 
to a shock, as the system now satisfies (B.ll I, 

B = D - y/\E(D)-E \\D-D \. 

For Do = 0, we have a shock in both directions, 

B = D - sign(D) v /|D(D)-D ||D-D |. 

For Do > 0, again, we must consider the different positions of D relative to the other points. These cases 
appear in Figure 18 If D < D 2 < 0, there is the shock solution satisfying (B.ll), 

D = Do + v / |D(D)-D ||D-D |. 
For D 2 < D < 0, this becomes a compound wave, 



D = D - / * y/E'(s)ds + y/\E(D)-E it \\D-D i< \ 



'D 

D* is again the point on the curve whose tangent intercepts (D,D(D)). For < D < Do, this becomes a 
purely rarefactory wave, 



f D 

B = D - / y / E'(s)ds. 

JDo 



! D 
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-E(D,z) 




D Do 



-E(D,z). 



Rightward Shock 
Wave 




(c) < D < D 2 



D Q D 



Rightward Rarefaclion 
Wave 



(b) D < D < 

-E(D,z)>- 



Rightward Shock 
Wave 




(d) D 2 < D 



FIGURE 17. Entropy satisfying rightward traveling waves when D < 0. 

Finally, for D > D , we again have a shock, 

B = B - y/\E(D)-E \ \D-D \. 
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(b) D 2 < D < 
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Rightward Shock 
Wave 
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(d) D < D 



FIGURE 18. Entropy satisfying rightward traveling waves when D > 0. 
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